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ABSTRACT 

Using  a  one  level,  barotropic  ocean  model,  driven  by  sur- 
face winds,  a  finite  difference  form  of  the  vorticity  equation 
was  integrated  over  210  days  of  simulated  time.   The  solutions 
using  constant  coefficients  of  lateral  eddy  viscosity  were 
compared  with  those  using  variable  coefficients  derived  from 
enstrophy  cascade  and  energy  cascade.   Using  a  constant  eddy 
viscosity  coefficient  of  rather  low  magnitude  produces  a  large 
amplitude  computational  oscillation  which  fills  the  entire 
basin.   An  order  of  magnitude  larger  coefficient  produces  a 
marginally  satisfactory  solution,  where  the  western  boundary 
current  was  rather  well  represented,  but  a  moderate  computa- 
tional oscillation  was  still  evident.   By  increasing  the  co- 
efficient yet  another  order  of  magnitude,  the  computational 
oscillation  is  negligible,  but  the  solution  in  the  ocean  in- 
terior is  unrealist ically  damped.   An  accurate  physical  and 
numerical  depiction  of  both  the  ocean  interior  and  western 
boundary  with  no  computational  oscillation  was  achieved  by 
using  either  of  the  two  forms  of  non  linear  eddy  viscosity. 
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I.   INTRODUCTION 

This  study  is  part  of  a  broader  effort  to  develop  the 
capability  of  making  large  scale  oceanographic  predictions 
on  a  dynamical  basis.   The  large-scale  and  meso-scale  thermal 
structure  of  the  ocean  is  very  dynamic  and  is  related  to  most 
marine  and  atmospheric  processes.   Some  of  the  most  obvious 
relationships  to  sea  surface  temperature  (SST)  are  ocean 
fronts,  circulation,  currents,  and  sea  state;  and  atmospheric 
temperature,  circulation,  and  wind  velocity.   The  effect  that 
these  environmental  phenomena  can  have  on  naval  and  maritime 
operations  is  well  known. 

The  significance  of  the  long  range  effect  of  SST  anomalies 
on  weather  patterns  has  received  increasing  scientific  aware- 
ness.  On  20  August  1974,  J.  Namias,  at  a  NORPAX  Conference 
at  the  U.S.  Naval  Postgraduate  School,  reported  on  the  results 
of  an  empirical  ocean/atmosphere  prediction  model.   With  a 
knowledge  of  SST  anomalies  in  the  spring  of  1974,  Namias 
derived  fields  of  atmospheric  temperature  anomalies  and  cir- 
culation for  the  following  summer.   Namias  predicted  in  May 
that  a  comparatively  severe  drought  would  occur  over  the 
midwestern  United  States  during  the  summer  of  1974.   Later 
events  verified  his  forecast. 

One  of  the  areas  of  difficulty  in  adequately  representing 
the  circulation  and  anomalies  in  a  finite-grid  ocean  circula- 
tion model,  and  subsequently  satisfying  and  integrating  the 


equations  of  motion,  is  the  representation  of  internal 
frictional  effects  within  the  ocean  fluid. 

Using  linear  theory,  Takano  [1974]  showed  that  if  the 
coefficient  of  eddy  viscosity  is  too  small,  a  computational 
oscillation  in  space  results  from  not  resolving  the  western 
boundary  current.   This  computational  oscillation  fills  the 
entire  basin  and  contaminates  the  solution.   On  the  other 
hand,  if  the  coefficient  is  increased,  the  open  ocean  solu- 
tion is  too  viscous.   To  handle  this  problem,  Takano  showed 
that  the  use  of  upstream  differencing  in  the  beta  term  of 
the  vorticity  equation  allows  a  smaller  coefficient  than  per- 
mitted by  linear  theory.   However,  this  method  unfortunately 
produces  excessive  damping  of  time  dependent  motion.   In 
addition,  this  mathematical  scheme  is  not  representative  of 
any  physical  motion  in  the  ocean.   A  better  approach,  from 
a  physical  point  of  view,  is  the  use  of  non  linear  eddy  vis- 
cosity.  The  friction  force  which  arises  using  non  linear 
eddy  viscosity  is  extremely  sensitive  to  the  scale  of  the 
motion,  and  therefore  will  be  relatively  small  in  the  oceans' 
interior  where  the  scales  are  comparatively  large,  and  will 
be  relatively  large  in  the  western  boundary  region  where  the 
scales  are  comparatively  small.   The  comparatively  large 
dissipation  in  the  western  boundary  region  will  keep  the 
current  broad  enough  to  be  resolved  by  the  grid  and  thereby 
prevent  the  formation  of  any  computational  oscillation. 

In  1968,  Leith  examined  two  dimensional  turbulence 
advection  and  derived  non  linear  coefficients  of  eddy 
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viscosity  from  the  cascade  of  enstrophy  from  large  scale  to 
small  scale  motion.   In  this  case,  the  coefficient  of  eddy 
viscosity  is  proportional  to  the  magnitude  of  the  horizontal 
gradient  of  vorticity. 

Another  satisfactory  non  linear  procedure  was  introduced 
by  Smagorinsky  [1963]  in  which  an  estimate  of  the  energy  cas- 
cade rate  is  made  from  the  fluid  deformation  itself.   In  this 
case,  the  coefficient  of  eddy  viscosity  is  proportional  to 
the  absolute  value  of  the  deformation. 

The  purpose  of  this  thesis  is  to  introduce  non  linear 
coefficients  of  lateral  eddy  viscosity  that  are  not  only 
reasonably  representative  of  actual  physical  conditions,  but 
these  same  coefficients  must  retain  their  reasonable  repre- 
sentation when  used  in  finite  grid  numerical  ocean  circulation 
models.   In  addition,  with  these  coefficients  the  numerical 
model  must  remain  computationally  stable  in  long  term  inte- 
gration.  Integrations  using  both  of  the  above  forms  of  non 
linear  eddy  viscosity  are  examined  and  compared  with  numerical 
and  analytic  results  using  constant  coefficients. 
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II.   THE  MATHEMATICAL  STATEMENT  OF  THE  PROBLEM 

A.   FORM  OF  THE  VORTICITY  EQUATION 

The  continuity  equation  for  a  homogeneous  fluid  and  the 
non  linear  equations  of  horizontal  motion,  cross  differenti- 
ated to  eliminate  the  pressure  term,  form  the  basis  for  the 
vorticity  equation: 


T„      3  T 


X 


dr  9       V     d 

d%  +  V8  +  fV  V  =   ^(Fy  +  p-^r)  -  gy(Fx  +  p^H) 

or      i^  +   V-Vc+v&    +    fv-V   =   AFy    +    T7   )   "  ^|(Fx  +  T*   ) 
9t  9*  PoH        dy  ^ToH 

(II-D 

j  r 

where  $  =  -5—  is  taken  constant.   In  (II  - 1)  ,  (Tx,  ?y)  is  the 
surface  stress,  H  is  the  constant  basin  depth,  and  all  other 
symbols  have  their  usual  meaning.   Bottom  stress  was  neglected, 

The  friction  forces  are  represented  by 

Fx  =  v(Avu) 

Fy  =  V-CAv  v)  (H-2) 

where  A  is  the  coefficient  of  lateral  eddy  viscosity  and  (u,v) 
is  the  vertically  averaged  velocity,  as  discussed  below. 

The  model  is  designed  with  a  "rigid  surface,"  in  order 
to  filter  gravity  waves  from  the  system: 

w(o)  =  0. 
It  is  also  necessary  that  there  be  no  vertical  motion  on  the 
flat  ocean  floor: 

w(-H)  =  0. 
Because  w  is  zero  at  the  top  and  bottom,  the  divergence  of 
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the  vertically  averaged  current,  V,  is  zero. 

V  •  V  =  0. 
Therefore  V  can  be  defined  by  a  streamf  unction ,  \\i   , 

u  =  -  |^ 
dy 

v  =   M 
8X. 

Integrating  (II-l)  from  top  to  bottom  and  using  y  V  =  0, 

the  final  form  of  the  vertically  integrated  vorticity  equation 

was  written 


dt    r      8x  .    3x  /    p0H/ 


x 


-  ^C-V-  VU  +  Fx  +  ^)  (n_3) 

Equation  (II-3)  also  includes  the  assumption  that 


V  • Vv  = V  •  v,  etc. 

The  next  section  will  discuss  the  formulation  of  the 
friction  terms,  Fx  and  Fy,  which  depend  upon  the  non  linear 
eddy  viscosity  coefficient,  A. 

B..  ENSTROPHY  CASCADE  FORM  OF  NON  LINEAR  EDDY  VISCOSITY 

The  first  method  of  generating  non  linear  coefficients 
of  eddy  viscosity  is  through  the  theory  of  two  dimensional 
turbulence,  in  which  enstrophy  and  kinetic  energy  are  con- 
served by  the  advective  terms.   Using  these  principles, 
Kraichman  [1967]  and  Leith  [1968]  derived  the  relation 

E(k)  oc  n2/3k"3  (II-4) 

where  r\    is  the  assumed  constant  cascading  rate  of  mean  squared 
vorticity  and  E (k)  is  the  kinetic  energy  in  wave  number  k. 
The  eddy  viscosity  which  causes  the  dissipation  is  also  assumed 
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a  function  of  n  and  k,  and  by  dimensional  analysis 

A  =  a  nV3k-2  (H-5) 

where  a  is  a  constant. 

One  estimate  of  n  made  locally  as  a  function  of  sur- 
rounding data  was  made  by  Leith 

n  =  AC  V  O  •  (  V  O  =  A|VC|2,  (II-6) 

where  £  is  the  local  vertical  component  of  vorticity. 
Substituting  (II -6)  into  (II-5)  and  solving  for  A  leads  to 

A  =  (  a  V*/k)  3|Vc|. 
Assuming  that  the  wave  number  k*  =  2iT/d,  where  d  is  the  grid 
size,  lies  in  the  inertial  subrange,  the  final  form  of  non 
linear  eddy  viscosity  for  enstrophy  cascade  becomes 

A  =  CaV*  _<L)  3|VC|  (II_7) 

Denoting  the  minimum  viscosity  coefficient  by  A0,  the 
viscosity  in  the  model  was  written, 

A  =  A0  +  Y|v^|d3  (11-8) 

A  is  a  maximum  when  |  V£  |  is  maximum.   Defining  Amax  =  mA0, 
then  (II-8)  becomes 

Amax  =  mA0  =  A0  +y|  Vc  Uax  d'  (H-9) 

Solving  for  y 

1  ^'max  a  (11-10) 


Y  = 


The   quantity    |Vc|max  was    estimated  by 

I^Uax   "    6.67xl0"llt         cm"1    sec"1  ("-ID 

which    for   d   =    300    km   leads    to    a  value    for  y    of 
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A0(m-1)  (II.1 

1.8xl09  J 

Substituting  this  into  (II -8)  gives  the  final  equation  for 

the  enstrophy  cascade  form  of  non  linear  eddy  viscosity  for 

this  study: 

A  =  A0(l  +  Chi-1)[Vc|)  =  ao(1  +  Qn-1)  |  Vg  I  ^  (II -13) 

I v? I  max  6-67 x  10 

In  (11-13),  Vmax   was  taken  from  linear  theory  [Munk,  1950], 
and  m  was  considered  an  adjustable  parameter. 

With  this  non  linear  form  of  eddy  viscosity,  the  vertically 
averaged  friction  force  terms  in  the  x  and  y  directions  become 
from  (II  -  2) 

Fx  =  -J-(A^)  +  _lCA2H)  (11-14) 

8x    3x     dy        dy 

Fy  =  MA^   +  MA^   •  (n-15) 

3         dx        9x     3y    dy 

The  vertical  component  of  the  curl  of  the  friction  force 

becomes 

CURL7  F  =i5x-8-I^  .  (11-16) 

9x    dy 

C.   KINETIC  ENERGY  CASCADE  FORM  OF  NON  LINEAR  EDDY  VISCOSITY 

From  diminsional  arguments,  Leith  [1968]  showed  that  in 
three  dimensional  turbulence,  if  the  energy  cascade  rate  from 
large  scale  to  small  scale  is  proportional  only  to  kinetic 
energy  dissipation  (e)  and  wave  number  (k) ,  then 

E(k)  =  ctTe^k"5/3.  (H-17) 

If,  in  addition  to  this  equation,  it  is  assumed  that  eddy 
viscosity,  which  produces  the  dissipation,  is  also  a  function 

of  e  and  k  only,  then  by  dimensional  arguments 
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A  =  a2  eV3k"V3  (11-18) 

where  k  lies  in  an  inertial  subrange.   Assuming  the  dissipa- 
tion rate  is  constant,  then  k  ~  V3  %  d 4/3 .   This  quasi-linear 
eddy  viscosity  is  dependent  only  on  grid  size,  and  hence 
latitude  if  the  grid  size  is  latitude  dependent.   But  more 
realistically,  dissipation  is  not  constant  nor  independent 
of  motion  or  grid  size.   Smagorinsky  [1963]  assumes  a  local 
value  of  dissipation: 

e  =  A| D | 2  (11-19) 

where  | D|  is  the  magnitude  of  the  deformation  tensor.   In 
this  case  (11-18)  becomes 

A  =  a23/2  |D|d2  (11-20) 


In  the  deformation  tensor  |D|  =  /  Ds2  +  D^2  ,  the  shearing 
deformation  is  Ds  =  QX.   +..JLIL.   and  stretching  deformation  is 


Dt  =  3'u  _  dv 

8x   37 


cTx   3y 


The  Smagorinsky  equation  takes  the  form  similar  to 
(II-8)  for  this  study 

A  =  A0  +  y  |  D|  d2  .  (H-21) 


Defining  Amax  =  mA0,  equation  (11-21)  becomes 

mA0  =  A0  +  Y|D|max  d2 

Solving  for  y 

Ao(m-l) 


Y 


Dlmax  d2 


(11-22) 


(11-23) 


Using  (11-23)  in  11-21),  the  final  energy  cascade  form  of 
the  non  linear  eddy  viscosity  becomes 
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=  An  [1  +  ISLiil 


D 
D 


] 


(11-24) 


max 


The  form  of  the  friction  force  principally  used  in  the 
model  for  the  kinetic  energy  cascade  case  was  of  the  form 
of  (11-14)  and  (11-15).   Experiments  following  the  friction 
force  of  Smagorinsky  [1963]  were  also  used  for  comparison: 


Fx  =JL(A  Dt)  +  -2-  (A  Ds) 
8x         dy 


(11-25) 

(11-26) 
In  both  cases  the  curl  of  the  friction  force  term  was  (11-16) 


Fv  =  JL.  (A  Ds)  -  -2-  (A  DO 

7      ax        zy  z 
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III.   THE  MODEL  AND  FINITE  DIFFERENCE  EQUATIONS 

A.   PHYSICAL  CHARACTERISTICS  OF  THE  MODEL 

The  model  consisted  of  a  one-level  barotropic  ocean  in  a 
square  basin  of  length,  L  =  9600  km;  of  breadth,  B  =  9600  km; 
and  a  flat  bottom  of  constant  depth,  H  =  2  km.   A  portion  of 
the  staggered  grid  is  shown  in  Figure  1.   The  total  grid 
points  (excluding  the  boundary  buffer)  are  33  x  33,  and  the 
distance  between  adjacent  grid  points  is 

X  =  Y  =  d  =  300  km.  (III-l) 

The  value  chosen  for  3  corresponds  to  a  grid  centered  at 
32.5°  N. 

The  staggered  grid  consists  primarily  of:  1)  the  inter- 
sections, or  principal  grid  points  (  •  ),  where  are  defined 
the  streamfunction  (  y  )  ,  vorticity  (  Z>  )  »  deformation  (  D  )  , 
and  the  coefficient  of  eddy  viscosity  (An)  generated  from 
deformation;  and  2)  the  grid  centers  (o),  where  are  defined 
velocity  (u,v),  the  gradient  of  vorticity  (v.c),  the  friction 
force  (Fx,  Fy) ,  and  the  coefficient  of  eddy  viscosity  (A^  ) 
generated  from  (Vc).   Auxiliary  variables  were  defined  at 
cross  (  x )  points . 

The  model  is  not  strictly  designed  to  simulate  any  partic- 
ular ocean,  but  rather  to  be  representative  of  the  fundamental 
physical  characteristics  of  mass  transport  and  western  boundary 
current  of  an  ocean  in  the  northern  hemisphere  as  affected  by 
the  viscous  action  of  the  ocean's  large  scale  circulation. 
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o 
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j+l    i+l,J+l     i+1, 
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-x- 
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j        i+l,j        i+1, 


j        i+l,j        i+1, 


1,1  1,1  2,1  2,\ 


Figure  1:     The  "tAggered  Grid  with  Corner  Boundary.     The  solid  line  represents 
the  physical  boundary  of  the  basin. 


(•  )     Principal  grid  point    (  \J>  ,    f  ,   D,   A(  |d|  )      ) 

(  o)     Grid  center   (u,  v,   V  s"    ,   A(  V  t,  ),  F   ,  F     ) 

(x)     Grid  average  points    (UAV,   VAV,    UAVE,  VAVE,   A(|d|  )ave,   T%  ,   Fyave    ) 
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The  energy  for  this  barotropic  wind  driven  circulation 
model  is  represented  by 

—  =  "  FcosCV-^'   Ty  ■  0,  (IH-2) 

a  pattern  of  westerly  winds  in  the  center  half  and  easterly 
winds  in  the  northern  and  southern  quarters  of  the  ocean 
(See  Figure  2) .   This  leads  to  a  stress  curl  term 

-2-fI*_<i   _-2nf   .  r2n^j,  Ax2  (IH-3) 

which  provides  the  actual  forcing  in  the  vorticity  equation. 
F  is  the  amplitude  of  the  zonal  component  of  the  wind,  and 
B  is  the  north-south  extent  of  the  domain. 

The  velocity  (u,v)  written  in  numerical  form  is 

"i,j  =  j[C^i-iJ-i  +  yj,j-i)  -  (!iiLillii)] 
vi,j  =     t  ^lliljLll^lzl        (^i-l,j  t  yj-i,j-D] 

2  2 

The  boundary  conditions  for  (u,v)  used  along  the  coastline 
were  zero  normal  flow  and  zero  slip.   Zero  normal  flow  was 
attained  by  requiring  the  streamfunction  ( ¥ )  on  the  left  side 
of  equation  (II-3)  to  be  equal  zero  on  the  ocean  perimeter. 
To  implement  the  condition  of  zero  normal  flow  and  zero  slip 
in  the  terms  on  the  right  hand  side  of  (II-3),  the  velocity 
is  defined  as  zero  on  the  coastline  by  defining  the  zonal 
boundaries 

Ui,l  =  -ui>2  ui,34  =  -ui,33 

Vi,l  =  -ui,2  vi,34  =  -vi,33,  (HI-5) 

for  the  meridional  boundaries 
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Figure  2:  Wind  Stress  Pattern 
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ul,j  ="u2,;  u34,j  ■  "u33,j 

Vl,j  ="V2,j  V34,j  "  -v33,j,  C111'6) 

and  for  the  corners 

Cu,v)1>]L  =  Cu,v)2>2     Cu,v)34>1  =  Cu,v)33>2 

Cu,v)1)34  =  (u,v)2j33    Cu,v)34)34  =  Cu,v)33>33    (III-7) 

where  the  original  32  x  32  (u,v)  grid  field  becomes  a 
34  x  34  grid  to  include  the  boundary  conditions. 

B.   FINITE  DIFFERENCE  FORM  OF  THE  VORTICITY  EQUATION 
The  finite  difference  form  of  (II-3)  was  written 

3t     i,J   P        2Ax 

i   (v-w)i+1  1  +  1  +  Ow)i+i,j 
-  TtC 2 ) 

(  V-Vv)itj  +  i  +  (v'Vv)i,j)] 

+  1  rf^'^iJ^  +  lV*Vu)i+l,j+K 
d  U         2  j 

(y-TaUi    +  (V.yu)i  +  1  j 

-   ( ^L_ LJ_  )] 

+  I  [(Fyj+i,j+i  +  Fn+i,j}  _  cFyj,j+i  +  Fyi>j)j 

Li  L\ 

-  1    r  fFxi,  i  +  1    +    Fxi-H.i  +  1-)     _     fFxi.i    +   Fxj+l,j 

3     LI  2  J  *.  2  J  J 

-  ^  sin    (^  yj)  d  (III-8) 
In    (III-8)    V2^.           is    given  by 


J 


'2*i,j  -jt  [»i*ifj  ♦Vi-u  ♦*i,j+i  ♦*i,,.i  -  4*i.i]- 
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The  advection  terms  (  V-Vu,  etc.)  were  finite  differenced 
using  the  same  method  as  Haney  [1974],  and  the  friction 
terms  were  expressed  differently  for  each  of  the  two  forms 
of  eddy  viscosity  as  discussed  below.   Readers  are  referred 
to  the  documented  computer  program  in  Appendix  A  for  further 
details . 

In  order  to  define  F,  it  is  first  necessary  to  determine 
three  values  of  minimum  A  (AMIN) ,  such  that  the  maximum  in 
the  analytic  streamfunction  (¥)  is  at  x  =  d/2,  d,  2d, 
respectively.   Since  the  grid  cannot  resolve  the  western 
boundary  in  the  first  case  (where  the  western  boundary  width 
is  d) ,  a  large  computational  mode  was  expected  with  AMIN(l). 
The  second  case  was  expected  to  be  marginally  stable  with 
AMIN (2);  and  finally,  with  AMIN(3)  it  was  expected  that  the 
western  boundary  would  be  clearly  resolved,  but  the  value  of 
AMIN(3)  would  be  so  high  as  to  unrealis tically  damp  the 
interior  ocean  solution.   The  respective  difference  equations 
for  the  three  AMINs  were: 

AMIN(l)  =  (3  x  {±     x  d/2) 3 

2n 

AMIN[2)  =  (3  x  ll  x  d) 3 

2n 

AMIN(3)  =  (3  x  {Z     x  2d)3  (III-9) 

2n 

These  three  values  of  AMIN  were  used  for  the  three  fundamental 
experiments  where  A  was  taken  as  constant,  and  also  for  the 
minimum  A  in  the  experiments  with  non  linear  coefficients  of 
eddy  viscosity  using  (11-13)  and  (11-24)  with  A0  =  AMIN. 
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In  determining  the  grid  size  and  location  for  the  non 
linear  coefficients  of  eddy  viscosity,  the  grid  location  of 
the  generating  parameters  (  |  Vc  |  or  |D|)  had  to  be  taken  into 
consideration.   This  led  to  a  34  x  34  grid  for  A(|Vc|)  in 
the  case  of  enstrophy  cascade,  and  a  33  x  33  grid  for  A(|D|) 
in  the  kinetic  energy  cascade  form.   (See  Figure  1)  . 

1 .   Enstrophy  Cascade  Case 

In  the  enstrophy  cascade  case,  in  order  to  generate 
coefficients  of  non  linear  eddy  viscosity,  it  is  necessary 
to  generate  a  relative  vorticity  field  and  to  determine  the 
gradient  of  vorticity  at  each  time  step.   Several  possible 
techniques  exist  in  developing  a  vorticity  field,  where  in 
all  cases  vorticity  would  be  defined  on  the  33  x  33  principal 
grid  points  (•)•   For  this  model  C  was  defined  in  terms  of 

(u,v)  by 

_  9v    8na 

ay 

VAV-  .  ,      VAV,      UAV,  ,  ,  -,    UAV. 


C    3x    ay 


^1,3  d  d  (111-10) 

i  =  1..  .33,  j  =  1..33 
where  UAV  and  VAV  are  defined  as 

Ui-i+Ui+li  Vii+Vi-i+l 

UAV-  •  -  U1»J  +  Ul  1»3   ;   VAV,  i  =  -^ ,J     , 

L»3  2  L»J         2 

(III-ll) 

i  =  1...33,  j  =  1...34;  i  =  1...34,  j  =  1..33 

in  order  to  get  an  equivalent  (u>v)  on  the  principal  grid 

points.   If  (u>v)  are  written  in  terms  of  Y,  then  (111-10) 

reduces  to 

Ci.j  -  1  Cn+lJ  +  1  +  *i+l,j-l  +  ^-1,3-1 

Zd  (111-12) 

H'i-l,j  +  l  -  4H<i  j) 
>j         »j  24 


This  method,  along  with  three  other  schemes  of  Miyakoda  [1962] 
for  calculating  vorticity,  were  utilized  and  compared. 

The  gradient  of  vorticity  C|Vc|)  was  developed  using 
centered  differences  on  the  32  x  32  grid  centers  ( ° } : 

lvUi,j    =    {[jUjiJ    +   gi,j-l    "   g  i  - 1 ,  .i    +   ^  i  - 1 ,  j  - 1 )  ]  ^ 

d  2  2 

+  [|(£LJ  +  EliLJ  "  Ci>J  +  ci-i,j-i)]2}^ 

2  2 

(111-13) 
i  =  2  ...  33  ,  j  =  2 , . . 33 

Finally,  (11-13)  in  finite  difference  form 
A(|Vc|)i  \    =  AMIN  x  (1  +  m  x  [V.gL-  a 

lV?Uax  (III  - 14) 

i  =  2. . . 33,  j  =  2. . .33 
where  Aj_  j  is  defined  on  the  32  x  32  grid  centers.   The  value 
of  [vcjm„   of  6  x  10 1 h   was  estimated  from  linear  theory 
[Munk,  1950],  and  m  was  considered  an  adjustable  parameter. 

The  finite  difference  form  of  the  friction  force  for 
enstrophy  cascade,  on  a  32  x  32  grid,  was  written  from  (11-14) 
and  (II-15)  as 

Fxi  i  =  1  r(Ai?j  +Ai  +  i,j)CUi+i,j  -  Ui,j)  - 

-  CAi-l,_i  +  Ai,j)(ui,j  "  Ui"l,j)] 

2 

+  1q  [(Ai,3  +  AiJ  +  i}{UiJ  +  i  -  Ui,j) 

-  fAi,J-l  +  Ai,J)CUi,J  "  Ui,J"D] 

2 

1  =  2. . .33,  j  =  2. . .33 
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Pat       -   1    Ai   -j   +  Ai+1   i 

Fyi,j  "  cPK   ?J  2     ,J)Cyi+l,J  -  vi,j) 

-  CAi-l,j  +  Ai,j)CVi,j  -  vi-ij)] 

+  ^[(Ai,j  +  Ai,j  +  luvi,3  +  l  -  vi,j) 

-  CAi,j-l  +  Ai,j)(vi,j  -  Vi,j-D]        (111-15) 

2 

v 

i  =  2. . .33,  j  =  2. . .33 
It  is  clear  that  (111-15)  requires  a  34  x  34  field 
of  A(|vc|)  in  order  to  calculate  the  friction  forces.   Two 
methods  were  utilized  to  generate  a  value  of  |Vc|  and  thereby 
a  value  of  eddy  viscosity  coefficients  across  the  boundary. 
The  first  case  was  linear  extrapolation  in  all  four  direc- 
tions.  For  the  western  boundary  the  finite  difference  form 

of  AC | Vc | )  was 

SLOPE  =  (A?  .  -  a,  ,)/d 

L  •>  J  °  >  J 

A-,  i   =  SLOPE  x  d  +  A,  .:  (III- 16) 

j  =  2 33 

and  the  values  of  A  immediately  outside  the  other  boundaries 
were  determined  in  a  similar  manner.   This  method  is  physically 
representative  because  it  gives  a  boundary  of  increasing  co- 
efficients in  the  direction  of  the  western  boundary.   Along 
the  other  boundaries  the  gradient,  and  hence  the  boundary 
value  of  A ( 1 V c | )  is  flat. 

The  second  method  utilized  to  generate  a  value  for 
A(|Vc|)  across  the  boundary  was  to  extend  the  existing  interior 
"boundary"  outwards  in  all  four  directions,  so  that  across  the 
western  boundary,  for  example 
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Al,j  ■  A2,j         3  =  2, ...33  (111-17) 

This  method  is  consistent  with  the  antisymmetric  velocity 
profile  which  accompanies  the  zero  flow  boundary  conditions. 
In  both  cases  the  34  x  34  corners  were  not  used  for  A(|Vc|). 
The  curl  of  the  friction  force,  defined  on  a  31  x  31 
interior  grid,  was 

CURL  Fi  i  =  CFyi'J  +  ¥-I±ilzl    -    Fyj-l>i  +  Fyj-i.1-1   ' 
i,j  2  2  Jta 

.  (Fxi-l,i  t    SiLJ  "  Fxi-lyj-l  +    FxiJ-l)/d 

i  =  3 33     j  =  3 33     (111-18) 

2 .   Kinetic  Energy  Cascade  Case 

The  first  steps  in  determining  the  deformation  field 
of  the  fluid  was   1)  to  calculate  the  shearing  deformation 

DSi  i  -  VAVj+lJ  -  VAVit1  +  UAVj.j.l  -  UAVj  ,  j 

>J  d  d  (111-19) 

i  =  1.  .33     j  =  1 33 

where  UAV  and  VAV  are  defined  in  (III-ll),  and  2)  to  calculate 

the  stretching  deformation  (Dt) 


Dt.  .  =  UAVEi  +  1J  -  UAVEifj    VAVEi>j  +  1  -  VAVEi>j 
1,3  d  d 

i  =  1.  .33     j  =  1. . . .33 
where  UAVE  and  VAVE  are  defined  as 

UAVE,  ,  =  ui,3  +  Ui,1  +  1  i   =    1---34 


(III  -20) 


i,J 


j  =  1 33 


VAVEi  .;  =  vi,i  +  Vj-H.i  i  =  1 33 

)J        ~2 j  =  1 34      (IH-21) 

Then  the  square  root  of  the  squares  of  the  values  of 
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deformation  along  and  normal  to  a  streamline  gave  the  value 
of  deformation  at  each  principal  grid  point  (•): 

l°li,j  =  [(Dsi,j)2  +  (^ij)2]1/2  (111-22) 

i  =  1 ....  33    j  =  1 , ...  33 

Next,  equation  (11-25)  in  finite  difference  form  was  used  to 

calculate  values  for  A(|d|)  based  on  kinetic  energy  cascade 

A(|D|)-  ,  =  AMIN(1  -  m  x  'Di>j  L)  (III-  23) 

1,3  |Dmax| 

i  =  1 ....  3  3   j  =  1,...33 
where  AMIN  and  m  are  as  described  above,  and  the  maximum 
deformation  (Dmax)  was  estimated  to  be  6.22  x  10"7.   Since 
A(|D|)  was  a  33  x  33  field,  the  form  of  the  friction  force 
did  not  require  additional  boundary  conditions  for  A(|D|). 

The  numerical  form  for  the  friction  forces  based  on 
A(|D|)  written  from  (11-14)  and  (11-15)  was 

„        1    Ai  -j  i  +  Ai  i 

F*i,j  =  5»[t  ,J    ,   '3)(U^l,3  -  Ui.j) 

z 

.  (Ai-l,j-l  +  Ai_iJ)(ui,j  -  ui-l,j)] 
2 

+  jpKAi,i  +  Ai-l,1)(ui,i^l  -  ui,j) 
2 

"  (Ai>3  +  ^±illl)(ui,j  "  u±, j-1)] 

2 

ifj  ■  y[(Ai»3-l  +  Ai,j)(vi^l,j  "  vi,j) 


Fy 


-  rAJ-i>j-i  *  Ai-i,juvi>j  "  Vi-1>J)] 

2 

+  ^2  [CAi,j  +  Ai-l,j)(viJ  +  l  -  vi,j) 
2 
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-    CAi,J-l    +   Ai-lJ-l)Cvi,i    -    vi,j-i)] 

CIH-24) 

The   Smagorinsky   form   of   the    friction    force    for   A  ^  |  D  |    from 
(11-26)    and    (11-27)    was 

Fxij^ECAi.j    ♦    Aifj-i)/2]     [CDtifj    +    DtiJ„1)/2] 

"     [(Ai-lJ    +   Ai_1J.1)/2]     [(Dti.lfj    +    Dti_1}j_1)/2] 
+     [(Ai,j    +    Ai-ifj)/2J     [(DSiJ    +'Dsi.ifj)/2] 
"     [CAi,j-l    +    Ai-i,j-i)/2]     [Dsi_1J.1)/2]} 

Fyifj-I<"Aifj  +  Ai,j_1)/2]  [(DSifj  *  DsiJ_1)/2] 
"  [CAi-ij  +  Ai_1J_1)/2]  [(Dsi_1J_1)/2] 
+    [CAij    +    Ai.1>j)/2]     [CDtifj    +    Dti.lfj)/2] 

"     [(AiJ-1    *    Ai-i,j-l)/2]     [D^.^    ♦    Dti.1J_1)/2]} 

( 1 1  I- 25) 
i   =    2. . . .33,      j    =    2 33 

The  finite  difference  form  of  the  curl  of  the  friction 

force  is  the  same  as  for  the  enstrophy  cascade  case  (111-18). 

3 .   Solution  Description 

Using  (III-8)  as  the  basis  for  solution,  a  centered 

time  difference  scheme  (Leapfrog)  was  used  for  all  terms 

except  the  friction  terms  which  were  evaluated  at  the  previous 

time  step.   With  time  steps  of  fourteen  hours,  equation  (III -8) 

was  integrated  for  210  days.   To  start  the  model  and  to  prevent 

solution  separation,  the  Euler-Backward  (Matsuno)  time  scheme 

was  utilized  for  the  first  and  every  fifty  time  steps.   In  the 
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kinetic  energy   cascade  experiments  the  form  of  the  friction 
force  shown  in  the  following  results  was  in  accordance  with 
equation  (111-24).   Experiments  were  also  conducted  using 
the  Smagorinsky  form  (III- 25)  which  gave  similar  results  to 
(111-24)  and  therefore  are  not  shown  in  this  study. 

The  solution  phase  of  the  model  had  to  solve  the 
equation 

V2(ff)  "  Fl  =  °  (111-26) 

for  the  tendency  3¥/9t.   This  was  done  using  a  direct 
Poisson  solver.   This  new  technique  was  written  by  R.  Sweet 
(1972],  based  on  a  method  originated  by  Buneman  [1969],  and 
revised  for  this  model  by  Professor  F.  Faulkner  of  the  Naval 
Postgraduate  School.   The  method  is  extremely  accurate  and 
made  possible  computer  time  savings  of  nearly  an  order  of 
magnitude . 
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IV.       RESULTS 

A.  RESTATEMENT  OF  PURPOSE 

The  main  purpose  of  this  thesis  was  to  present  a  scheme 
whereby  the  accuracy  of  the  numerical  solution  of  a  finite 
grid  ocean  circulation  model  would  be  improved  by  the  intro- 
duction of  non  linear  lateral  eddy  viscosity  coefficients. 
As  shown  in  earlier  chapters,  the  gradient  of  relative 
vorticity  or  the  fluid  deformation  could  be  used  as  the 
respective  parameters  to  generate  the  coefficients  of  eddy 
viscosity.   The  solutions  using  constant  coefficients  of 
lateral  eddy  viscosity  will  be  compared  with  those  using 
variable  coefficients  derived  from  enstrophy  cascade  (A^|vc|) 
and  kinetic  energy  cascade  (A ^  |D|)  respectively. 

B.  ANALYTICAL  CONSIDERATIONS 

The  initial  experiments  investigated  three  cases  of  con- 
stant coefficients  of  eddy  viscosity.   First,  an  accurate 
analytical  solution  for  the  streamfunction  was  made  by  means 
of  a  separate  model  where  grid  spacing  was  60  km.   The  ana- 
lytical solution  of  the  streamfunction  in  terms  of  the  wind 
stress  curl  was  developed  by  Munk  [1950],  who  considered  a 
linear  eddy  viscosity: 

Kx,y)  =  -rX(x)S_1CURLzT  (IV-1) 

where  r  =  domain  width,  curl  t  =  the  k  component  of  the  wind 
stress  curl,  and 
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X(x)  =  Ke-Y^x  cosC^I  kx  +  j£f  -  f-)    +  1 

-  -L  (kx-e"1^  "^  -  1) 
kr 

in  which  X(x)  =  distance  eastward  from  the  western  boundary, 
k  =  C3/A)  V3      and    K   =   -~z   -    ^ 


The  reader  is  referred  to  Figure  3  which  portrays  the 
analytical  s treamfunction  (  \p  )  for  the  three  cases  where  the 
maximum  of  \p   occurs  at  d/2,  d  and  2d,  respectively.   The 
three  values  of  A  which  permitted  the  above  analytical  situ- 
ation to  occur  were  the  three  constant  coefficients  of  eddy 
viscosity  now  examined  in  the  numerical  model. 

C.   EXPERIMENTS  WITH  CONSTANT  COEFFICIENTS  OF  EDDY  VISCOSITY 

The  first  constant  coefficient  of  eddy  viscosity, 
Aj  =  0.12  x  10 8  cm2  sec"1,  which  physically  represents  the 
interior  ocean  circulation  most  accurately,  produces  a  large 
amplitude  computational  oscillation  which  fills  the  entire 
basin  of  the  numerical  model.   Figure  4  shows  the  extent  of 
the  oscillation  produced  in  the  ip  field  by  this  relatively 
low  magnitude  coefficient  of  eddy  viscosity.   Figure  7  shows 
a  direct  comparison  of  the  analytical  ip  field  and  the  numer- 
ical \p   field  produced  by  Aj  at  the  latitude  of  maximum  wind 
stress  curl.   The  reader  is  also  referred  to  Table  I  which 
presents  a  tabular  comparison  of  \p   field  highlights  as 
generated  by  constant  eddy  viscosity  coefficients.   All  of 
these  features  are  in  accordance  with  the  study  by  Takano 
[1975]  . 
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Figure  4:   Numerical  Solution  of  Streamfunct ion 
for  A  !  =  Constant 
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Next,  A2  =  0.93  x  108  cm   sec"  ,  an  order  of  magnitude 
larger  coefficient  of  eddy  viscosity,  was  examined  in  the 
numerical  model.   Figure  5  shows  the  circulation  for  constant 
A2 ;  Figure  11,  a  graphical  comparison  of  analytical  and 
numerical  ty   field  for  A2 ;  and  Table  I,  a  tabular  comparison. 
A2 ,  used  as  a  constant,  produced  a  marginally  satisfactory 
solution,  where  the  western  boundary  current  was  rather  well 
represented,  but  a  moderate  computational  oscillation  was 
still  evident  with  the  value  of  the  maximum  streamfunction 
50%  higher  than  the  analytical  solution.' 

By  increasing  the  coefficient  yet  another  order  of 
magnitude  to  A3  =  7.5  x  108  cm2  sec"1  (Figure  6,  Figure  14, 
and  Table  I) ,  the  computational  oscillation  is  negligible 
with  only  about  3%  error  in  the  numerical  solution.   The 
western  boundary  current  placement  was  in  accord  with  the 
analytical  case  (Figure  3) ,  but  the  solution  in  the  ocean 
interior  was  unrealistically  damped  by  the  large  viscosity. 

D.   EXPERIMENTS  WITH  NON  LINEAR  COEFFICIENTS  OF  EDDY  VISCOSITY 

It  is  clear  that  there  is  no  single  coefficient  of  eddy 
viscosity  that  can  physically  or  numerically  depict  all  the 
aspects  of  fluid  circulation  both  in  the  ocean  interior  and 
in  the  western  boundary.   As  can  be  observed  from  the  results 
so  far,  the  objective  of  using  non  linear  eddy  viscosity  is 
to  have  low  coefficients  of  eddy  viscosity  in  the  ocean 
interior,  and  increasing  coefficients  approaching  the  western 
boundary  in  order  to  resolve  the  western  boundary  current  and 
to  prevent  the  development  of  a  computational  mode  in  the 
solution.  35 
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Figure  5:   Numerical  Solution  of  Streamfunct ion 
for  A2  =  Constant 
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Figure  6:   numerical  Solution  of  Stremfunct ion 
for  A 3  =  Constant 
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Non  linear  coefficients  that  represent  actual  physical 
processes  in  the  ocean  were  developed  in  Chapters  II  and 
III.   It  is  shown  below  that  if  the  limits  of  the  range  of 
these  coefficients  are  properly  chosen  the  desired  objective 
is  achieved. 

In  the  first  non  linear  experiments,  ranges  of  coeffi- 
cients of  eddy  viscosity  were  chosen  with  the  minimum 
coefficient  equal  to  AMIN(l)  =  A} ,  and  the  range  of  variation 
of  the  coefficients  was  governed  by  the  adjustable  parameter 
m  in  both  (11-13)  for  enstrophy  cascade  and  (11-25)  for 
kinetic  energy  cascade.   Of  course,  in  this  case  where  AMIN 
is  the  smallest,  the  greatest  range  of  m  was  required  to 
properly  resolve  the  western  boundary  and  to  prevent  an  un- 
acceptable computational  oscillation  from  developing  in  the 
solution.   Experiments  were  conducted  with  m  varying  from 
10  to  1000  with  results  that  are  noted  below.   It  should  be 
noted  that  the  value  for  m  is  not  precisely  a  direct  multi- 
plier for  the  range  of  A  due  to  the  non  linear  effect  of 
|Vcmax|  in  (IH-14)  and  |D|max  in  (II I  - 15)  .   The  actual 
values  of  A/AMIN  were  printed  out  in  the  experiments  and  the 
range  of  A/AMIN  appears  in  Table  II  for  the  most  significant 
experiments . 

In  the  enstrophy  cascade  (A  ^  |  Vc;  | )  experiment  for 
Aj  =  0.12  x  10 8  and  m  =  10  and  using  the  symmetrical  boundary 
conditions  for  A  as  indicated  in  (111-17),  the  actual  range 
of  A/Aj  was  1.0  to  7.0.   The  resulting  ty  field  can  be  seen 
in  Figure  8A.   Using  the  linear  extrapolation  shown  in 
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Figure  8A:   Numerical  Solutions  of  V  for  Ai  and  m 
Symmetric  Boundary  Conditions  for  A%|Vc 
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(111-16)  to  obtain  A  across  the  boundary,  the  solution  shown 
in  Figure  8B  was  obtained.   In  this  case  the  range  of  A/Ax 
was  1.2  to  7.9.   It  is  noted  here  that  in  the  extrapolation 
boundary  condition  the  maximum  A  occurs  outside  the  western 
boundary  and  therefore  is  artificially  derived;  and  that  the 
A/Ax  minimum  values  indicated  in  these  experiments  are  not 
actually  the  lowest  ratio  obtained,  but  a  value  representative 
of  the  A  obtained  in  the  interior  ocean.   In  the  kinetic  energy 
cascade  experiment  (A  a.  |d|)  for  this  case,  the  range  of  A/A 
was  1.1  to  10.8,  with  results  very  similar  to  the  enstrophy 
cascade  experiments  (Figure  8C) .   As  can  be  readily  seen  in 
the  figures  for  the  above  three  experiments  with  AMIN  =  k1 
and  m  =  10,  a  heavy  computational  oscillation  is  still  very 
much  in  existence,  although  a  definite  improvement  over 
Aj  =  constant  (Figure  4)  is  apparent. 

It  can  also  be  noted  here  and  in  the  following  experi- 
ments, that  although  the  method  used  in  deriving  the  non 
linear  coefficients  of  eddy  viscosity,  namely  A^  |Vc|  and 
A'v  |  D  |  ,  show  somewhat  different  characteristics  in  the  ocean 
circulation  pattern,  the  results  were  analogous  enough  that 
the  main  purpose  of  this  thesis  could  have  been  accomplished 
with  either  method  and  either  boundary  conditions  for  A  in 
the  case  of  A^|Vc|  .   In  addition,  several  tests  using  various 
combinations  of  the  Laplacian  given  by  (111-12)  and  the  usual 
5-point  Laplacian  were  made.   It  was  found  that  the  method  of 
defining  the  relative  vorticity  field  was  of  no  consequence 
in  the  results  achieved,  and  therefore  this  paper  does  not 
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Figure  8B:   numerical  Solutions  of  ¥  for  Ax  and  m  =  10 
Extrapolated  Boundary  Conditions  for  A^jVc 
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warrant  any  further  discussion  of  this  aspect.   In  all  the 
above  non  linear  experiments  with  A3 ,  in  addition  to  the 
figures  listed  and  Table  II,  the  reader  is  referred  to  the 
graphical  representation  of  the  steady  state  maximum  analyt- 
ical and  numerical  streamfunction  field  in  Figure  7. 

In  the  next  set  of  experiments  with  AMIN(l),  m  is  in- 
creased to  100  giving  a  two  order  of  magnitude  range  for  the 
coefficients  of  eddy  viscosity.   The  results  are  shown  in 
Figures  9A,  9B ,  9C,  7,  and  Table  II.   In  the  cases  where 
A^|Vc|  with  the  symmetric  boundary  conditions  for  A,  A/A, 
ranged  from  1.2  to  41.7,  and  for  the  linear  extrapolation 
boundary  conditions  A/A.j  ranged  from  1.3  to  62.2.   For  A^|D| 
the  values  of  the  coefficients  ranged  from  1.5  to  70.1.   In 
this  group  of  experiments  the  solution  of  streamfunction 
improved  greatly,  with  the  numerical  ^max  within  a  few  per 
cent  of  the  analytical  ^max  in  all  cases.   However,  the  higher 
value  of  A  at  the  western  boundary  resulted  in  a  tendency  for 
^m^Y  to  move  eastward,  especially  in  the  case  of  A^|d| . 

Increasing  the  range  of  Aj  another  order  of  magnitude, 

the  next  experiments  examined  m  =  1000.   These  results  are 

shown  in  Figures  10A,  10B,  10C,  7,  and  Table  II.   In  all  cases 

the  computational  oscillation  was  negligible,  the  interior 

solution  was  not  strongly  damped,  and  the  western  boundary 

region  was  well  defined.   The  variation  of  A/Aj  was  1.5  to 

165  for  A%|vc|  (symmetric  boundary  conditions),  1.6  to  237 

for  A^|Vc|  (extrapolated  boundary  condition),  and  5.0  to  313 

for  A^lDl.   Again,  however,  there  was  a  tendency  for  i(;_„v  to 

max 
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Figure  9A:   Numerical  Solutions  of  y  for  Ai  and  m  =  100 
Symmetric  Boundary  Conditions  for  A^|V^| 
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Figure  9B:   Numerical  Solutions  of  y  for  Aj  and  m  =  100 
Extrapolated  Boundary  Conditions  for  A%|vH 
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Figure  9C 


Numerical  Solutions  of  y  for  A!  and 
m  =  100,  A^|D| 
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10A 


Figure  10A:   Numerical  Solutions  of  y  for  Ax    and  m  =  1000 
Symmetric  Boundary  Conditions  for  A%|Vr,J 
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Figure  10B 


numerical  Solutions  of  y  for  Ai  and  m  =  1000 
Extrapolated  Boundary  Conditions  for  A^lvH 
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Figure  IOC 


Numerical  Solutions  of  y  for  Aj  and 
m  =  1000,  A%|D| 
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move  eastward  to  3d   in  the  case  of  A<\,  |  D I  .   Due  to  this,  and 
the  slightly  higher  minimum  coefficients  in  the  ocean  interior 
for  the  A'v  |  D  |  solution,  the  A*v|vc|  solution  appears  to  be 
generally  preferred  to  A^IDI.   This  is  probably  because  the 
method  based  on  two  dimensional  turbulence  is  somwhat  more 
sensitive  to  the  space  scale  of  motion. 

In  the  next  sets  of  experiments,  AMIN  was  increased  to 
A2  =  0.93  x  10 8  cm2  sec"1  and  m  was  examined  for  values  ranging 
from  10  to  100.   These  results  are  shown  in  Figures  11,  12A, 
12B,  13,  and  also  Table  III.   With  a  higher  initial  AMIN  and 
resulting  higher  viscous  solution  in  the  ocean  interior,  the 
computational  oscillation  was  suppressed  and  a  completely 
satisfactory  solution  was  attained  by  m  =  20.   The  range  of 
values  for  A/A  and  \bmct„   are  given  in  Table  III,  and  Figure  11 
gives  a  graphical  representation  of  results  for  the  non  linear 
experiments  with  A2 .   Increasing  m  to  100  had  the  undesirable 
effect  of  moving  the  western  boundary  eastward  to  3d   in  the 
case  of  A^|D | . 

Non  linear  experiments  with  AMIN  =  A3  =  7.5  x  10 8  cm2  sec"1 
produced  no  enhancing  results.   The  solution  was  exceptionally 
viscous,  and  the  interior  ocean  was  already  overdamped  with 
such  a  high  minimum  coefficient  of  eddy  viscosity.   Increasing 
m  to  10  resulted  in  moving  the  tyma„    for  A^  |  D|  eastward  to 

JU  cIa 

3d.    Refer  to  Table  III  for  representative  values  of  A/A3 
and  ^max,  and  to  Figure  14  for  the  streamfunction  field  for 
A  'v,  I  D I  with  m  =  10. 
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Figure  12A:   Numerical  Solutions  of  <F  for  A2  and  m  =  10 
Extrapolated  Boundary  Conditions  for  A 
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Figure    12B 


ilumerical    Solutions    of    H'    for   A2    and 
m    =    10,    A%|D| 
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Figure  13 


numerical  Solutions  of  ¥  for  A2  and  m  =  21, 
Extrapolated  Boundary  Conditions  for  A 
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The  author  is  of  the  opinion  that  the  experiment  with 
AMIN  =  Aj  (or  perhaps  a  slightly  higher  minimum  coefficient) 
and  m  ranging  somewhere  between  100  and  1000  (depending  on 
minimum  acceptable  computational  oscillation)  gives  the  best 
field  of  non  linear  lateral  eddy  viscosity  coefficients  for 
future  utilization  in  more  sophisticated  finite  difference 
ocean  circulation  models.   It  appears  from  the  experiments 
that  the  only  drawback  is  that  these  higher  values  of  m 
produce  a  western  boundary  width  which  is  nearly  2d. 

An  accurate  physical  and  numerical  depiction  of  both  the 
ocean  interior  and  the  western  boundary  with  no  computational 
oscillation  was  achieved  by  using  either  of  the  two  forms  of 
non  linear  eddy  viscosity.   These  were  achieved  with  the 
minimum  coefficients  approximately  equal  to  Ax    in  the  interior 
ocean,  increasing  approximately  two  orders  of  magnitude  in 
the  more  dynamic  flow  in  the  western  boundary  region. 
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Hi 


Figure  14:   Numerical  Solutions  of  f  for  A3  and 

m  =  10,  A%| D| 
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V.   CONCLUSIONS 

This  numerical  model  has  been  successful  in  testing  a 
scheme  whereby  the  accuracy  of  the  numerical  solution  of  a 
finite  grid  ocean  circulation  model  can  be  improved  by  the 
introduction  of  non  linear  lateral  eddy  viscosity  coefficients. 
Non  linear  coefficients,  properly  generated  to  represent 
actual  physical  processes  that  may  be  occurring  in  the  ocean, 
allow  the  use  of  low  coefficients  in  the  interior  ocean  solu- 
tion, and  high  coefficients  in  the  higher  circulation  density 
of  the  western  boundary  current.   This  distribution  of  eddy 
viscosity  is  sufficient  to  prevent  the  formation  of  a  compu- 
tational oscillation  which  would  occur  if  the  low  value  were 
used  throughout  the  domain. 

Two  methods,  namely  enstrophy  cascade  (A^  |VC|)  and  kinetic 
energy  cascade  (A ^  |D|),  were  investigated  and  presented  which 
will  allow  the  researcher  to  generate  non  linear  coefficients 
in  other  models.   From  experimentation  with  this  barotropic 
model,  minimum  coefficients  and  corresponding  ranges  of 
coefficient  magnitudes  are  recommended  for  initial  investiga- 
tion in  more  sophisticated  baroclinic  ocean  and  coupled 
ocean-atmosphere  finite  difference  models. 
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APPENDIX  A 
C 

c 

C  LCOR  J.  M.  WPIGHT,  JP  USN.   XM-34,  .  DEPT  CF  METEOPOLOoY 

C  U.  S.  NAVAL  POSTGRADUATE  SCHPOL,  MCNTFREY,  CALIFORNIA 

C  THESIS:   OCEAN  CIRCULATION  M3DEL,  ONE  LEVEL  BAPOTROPIC  PHASE 

C 

C 

C  THE    MAIN    PROGRAM     INITIATES    THE     SOLUTION    OF    THE     VO°TICITY 

C  EQIATT.CN    BY     ZEROING    ALL    DIMENSIONED    V AR I ABL E S . THE    S  DfiROUT  IN  E  S  t 

C  NAMELY    CALF    AN!)    SOLVE",     THEN    TAKE    CONTROL    AND    PERFORM    THEIR 

C  RESPECTIVE     FUNCT I ONS . A F^ E RWARD,     THE    MA  I  \    PROGRAM    TAKES    CONTROL 

C  AND    PERFORMS     THE     TIME-STEPPING    PHASE    OF    THE     SOLUTION. 

C 

C  CVAR      MEANS   DIMENSIONED  VARIABLES     *********************** 

C 


REALMS  TITLEK12) 

REAL  *4  CL(41) 

LOGICAL* 1  LTG( 3) /.TRUE. , .TRUE. i .FALSE./ 

I  M  =  3  3 

JM  =  33 

I  MM  1=32 

JMM1=32 

T1P1=34 

JMP1=34 

CL  (  U  =  -4.0 

DO  10  N=2,41 
10  CL(N)=CL(N-l)*-0.2 
C 

C      STAPT  THE  ZEROING  PROCESS. 
C 

03  3C00  J=1,JM 

00  3000  1=1, IM 

PS  [Mil  I  ,  J)=0.0 

PSI(I»J)=0.0 

F1(I,J)=0.0 

RESIDI I, J)=0.0 

PSTEMtM  I  ,J  )  =0.0 
3C00  DPS!CT( !  ,  J)=0.0 

00  4000  J=1,JMM1 

DO  4000  1  =  1,  IMV1 

DELSGV( I  ,  J)=Q.O 
4000  DELS0U( I  , J  )=0.0 

DO  5C00  J=1,JMP1 

00  5C0O  I=1,IMP1 

A{  I,  J)=0. 

U(  !,  J)  =  3.0 
5C0O  VI  I,  J)=).0 

KK  =  0 
C 

C    MATSUNU  SCHEME  IS  USED  FOR  THE  FIRST  TIME  STEP,  AND  EVERY  50  TIME 
C       STEPS  THtPEAFTER 
C 

MATSNO=5  3 

DELTAX=3.0  *  10.0**7 

DELTAT=  5.0  *10.0**4- 

TIMEST=210.C*B6400. ) 

TI^E=0.0 
C 
C      LFAPFPOG  TIME  SCHEME   ****************************************** 

C 
C 
C      THE  LEAPFROG  SCHEME  IS  NEUTRAL  IN  CHARACTER  BUT  THE  PRESENCE  AND 
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C  UTILIZATION  OF  THREE  TIME  LEVELS  IN  THE  DESCRIPTION  OF  THE  FIRST 

C  DERIVATIVE  C'F  t>S  I  WITH  RESPECT  TO  TIME  ,  PRODUCES  A  COMPUTATION- 

C  AL  MODE  IN  TIME. THIS  MODE  HAS  TO  BE  AND  IS  REMOVED  bY  THE  EULER- 

C  BACKWARD  SCHEME. THIS   SCHEME  USES  TwU  LEVELS  IN  TIME  IN  ITS  UE- 

C  SCRIPTIPN  CF  THL  FIRST  T  I  M  F  DERIVATIVE  OF  PSI.  THEF  E  FOR  E  ,  THE  P.  E 

C  IS  NO  C  JMPUTATIUNAL  MOOF  IN  TIME  AND  A  '^ORE  ACCURATE  SOLUTION 

C  FOR  THE  DSI  FIELD  AT  tjm.f  LEVEL  'N+l'  IS  ATTAINED  .   *********** 


C 


DC  3100  K=lr3C00 

IF(K  .E  3.  1  )   GO  TO   2000 

L=MC'C(K,  MATSrjQ) 

IF  (L.LO.O)   GO  TO  2000 

CALL  CALF 

CALL  SOLVER (DPS IDT, IMM1, JMMl , DELTAX, DELTAX) 

DO  6000  J=2,JM"1 

DO  6CCC  1=2, IMM1 

TEMP=PST ( I , J) 

PSI!I,J)=  PSIMK  I  ,  J) +2.0  *  DELTAT  *  DPSIDT(I,J) 
6000  PSIMHI.J)  =  TEMP 
C 

C      PRINT  THE  PSI  FIELD  IN  TABULAR  FORM  **************************** 
C 

GO  T0   2300 
C 
C 

C  AT  THIS  time  ,  psi  IS  KNOWN  AT  TWO  CONSEQUETIVE  TIME  LEVELS  AND 
C  IS  STOFJ  0  H  PS!  AND  PSIMl  FIELDS  RESPECTIVELY.  BY  MEANS  OF  THE 
C  TEMP  STATEMENT  ,  THE  TWO  TIME  LEVELS  ARE  ADVANCED  BY  ONE  IMCP.E- 
C  MENT  OF  TIME  ,  I.E.,  DELTAT  =  5.0  *10.0**4  <SEC.)  RESPECTIVELY. 
C 
C 
C 

C      THE  LULER-6ACKWARD  OR  MATSUNQ  TIME  SCHEME  ;  NOTE:  PSI  AT  TIME 
C      LEVEL  »N«  IS  SAVED  IN  »STEMP  FIELD  FOP  LATER  USE  IN  SUBSEQUENT 
C      BACKWARD  IMPLICIT   STE^1  .   ************-.************************ 
C 
C 
2000  WRITEI6,  103)  TIME 
2100  00  2200  J=2,JMM1 

DO    2200     1=2, INM1 

psinu  i ,  J)  =psi  (i ,  j) 

2200  PSTFMPl I , J)=PS!( I, J) 
C 
C 

C      BY  MEANS  OF  THE  ASbVE  LOOP ( *2200 ), PS  I  AT  TIME  LEVEL  N  IS  PUT  IN- 
C      TO  THRE'i  FIELDS  .NAMELY  PSI.PSTEMP  AND  PSIMl.  ALL  FIELDS  ARE  AT 
C      THE  SAMF  time  LEVEL  IN  3FDER  TO  PRESERVE  LINEAR  COMPUTATIONAL 
C      STABILITY  IN  THE  FRICTION  TERM  OF  THE  FORCING  FUNCTION, Fl.  ***** 
C 

C      CC^MENCE  TEE  FORWARD  TIME  STEP  PHASE  OF  THE  MATSUNO  SCHEME.  **** 
C 

C      NOTE:  ALL  TERMS  OF  CALF  AND  RELAX  APE  OUTPUTED  AT  TIME  LEVEL  N  . 
C 

CALL  CALF 

SCALE=10.0**(-2) 

DO  2910  J=2,JMM1 

JS=JMP1-J 

DO  2910  I=2,IMM1 
2910  RESID(  I, J)  =  F1(  I, JS)*SCALE 

WRITE(6,138)     <(PESIO(T,J),I=2,32),J=2,32l 
103    FORMAT! • C«  ,'RESID    FIELD    =     Fl    FIELD    *    S  CALE  ■  ,  /  ,  3  1  (  1  X,  3  1F<V 


1.2, //I  ) 

CALL    SOLVER  I  DPS  IDT, I  MM  1 , JMMl , DELTAX, DELTAX ) 


C 

C      DDSIDT  IS  NOW  AT  TIME  LEVEL  'N'  .   ***************************** 

C 

DO  2500  J=2,JMMi 

DO  2500  1=2, IMM1 

PSI(I,J)=  PSIII.J)  ♦  DELTAT  *  DPSI0T(I,J) 
2500  PSIMl! I , J)=PSI  II , J) 
C 

C      PRESENTLY,  THE  FIELDS  PSI, PSIMl  AND  PSTEMP  HAVE  CONTAINED  IN   ** 
C      THEM  ,  PSI  VALUES  AT  TIME  LEVELS  'N-H'  INTERMEDIATE  ,'N+l'INTER- 
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MEDIATE    AND    *N*     RESPECTIVELY.        ********************************** 


700 
800 


NGU  AS  A  PART  OF  THE  MATSUNO  TI^E  SCHEMF  ;  A 
E)  BACKWARD  TIME  STEP  IS  EMPLOYED.  THE  BACKW 
A  CALIBRATION  STEP  TO  REFINE  THE  PSI  FIELD  A 
COMMENCE  TEE  BACKWARD  TIME  STEP  PHASE  OF  THE 
NOTE:  ALL  TERMS  OF  CALF  CALCULATED  AS  A  FU\'C 
MEDIATE  TIME  LEVEL  '  N+- 1 '  .   ***************** 

CALL  CALF 

CALL  SOLVER(DPSIDT, IMM1, JMM1,DELTAX,DELTAX) 

CPSICT  AT  TIME  LEVEL  •  N*  1  •  (INTERMEDIATE)  "U 
NOV,  OUTPUT  PSI  AT  TIME  LEVEL  'N+l*.  ******** 

DO  2700  J=2,JMM1 
DO  2  700  I=2tIMMl 
PSIMK  I,  J)=PSTE*P(  I  ,  J  I 
PSI ( I, J)=PSTEMP( I, J)  + 
TIME  =  DELTAT*K 
WRITE(6,103)  TIME,K 
SCALE=  10.0**(-8) 


SIMULATED  (IMPLICIT 
ARD  STEP  IS  USED  AS 
T  TIME  LEVEL  •  N  ♦- 1  •  . 

MATSUNO  SCHEME  .  ** 
TION  OF  PSI  AT  INTFR- 
******************** 


TPUTED.   *********** 
******************** 


DELTAT  *  DPS  I DT ( I  , J ) 


;      PRINT  OUT  «A.«  FIELD  DERIVED  FROM  OFFORMATION 

DO  1140  J=1,JMP1 
JS=JM+2-J 
DO  1  140  1  =  1,  IMP1 
1140  RESID(  I, J)=A(  1  ,JS)/AMIN 
IF  (M0D(K,20) .EU.OJ 
aWRITE  (6  ,1145)  ((RESI C(  I, J  )  ,1  =  1  ,33) , J  =  l ,34) 

1145  FORMAT ( '0' ,' RESID  FIELD  =  A  FI ELD/AM  IN ',/, 34 ( 1 X, 33F4 . 0, //) ) 

DO  2900  J=1,JM 

JS=J^P1-J 

DO  29C0  1  =  1,  IM 
2900  RESIC(I,J)=  PSKI.JS)  *  SCALE 

IF     (MOD( K, 20)  .EG.O) 
3WRITE(6,  104)        ((RESID(I,J) 


103 
04 


1=1,33) ,J=1,33) 
FOPMATP     «,Tfc3,'TIME  =  ',F10.1,l6,//) 
FOPMATCO'  , 'RESID    FIELD    =    PSI     FIELD    *    SCALE 
LL=MCD(K,3) 
TDAY=TIME/£64C0.0 
WPITF(6,103)     TCAY.LL 

LCCP     IS    FOR    THE    CONVERSION    OF    DATA    TO    A    FORM 
-TUR'       .**************************<•********* 


2950 


>970 
1999 

980 
UOO 


,/,33(lX,33F4.2,//) ) 


MORE    USEABLE     IN     'CC 
********** 


DC    2950    J=1,JM 

DO    2950     1=1,  IM 

DATAH  I  ,  J)  =PSI  (I  ,  J)*SCALE 

IF ( (TDAY.GT. 202.0)  .AND. (LL.EQ.O)  )  GO    TO    9970 

GO    TO       9980 

RE/SL  (5,9999)     (  T  I  TLE  1  (  J  )  ,  J  =  l ,  12  ) 

F0PMAT(6A8) 

CALL    CCrjTUR(CATAl,33,33,33,CL,41,TITLEl,8,08,LTG) 

IF  (TIME.GE.TIMEST)     STOP 

CONT  INUE 

STOP 

DEBUG 

AT    3100 

TRACE    CN 

DISPLAY    K 

END 
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SUBROUTINE  CALF 

THIS  SUBROUTINE  CALCULATES  THE  RIGHT  HAND  SIDE  OF  THE  VORTICITY 
EQUATION  WHICH  IS  MULTIPLIED  BY  thf  SQUARE  OF  DELTAX. 


DVAR 


MEANS   DIMENSIONED  VARIABLES     *********************** 


CCMMCN/OVAR/PSIM1.PSI  ,  Fl  ,U,V,DELSQU,  DELS'JV,DPSIDT,RESID,PSTEMP,A 
COMMON/DOMAIN/ IM, JM, IMM1 . J^Ml , IMPl ,  JM°1 , K, KK.AMI N,DFF 
DIMENSION  PSIM1(33,33),PSI(33,33),F1!33,33),U(34,34),V(34,34), 
1DELSCU(3  2,  32)  .OELSQV  (32,  32)  ,  DPSIOTJ  33,  331  ,P.ESID(  34,  34)  , 
2PSTEKP< 33, 33)  ,PS  IANL <  100) ,  DEF(33,33t,  CEFS(33,33),  DEFT(33,33) 
DIMENSION  OELVGP (33, 33),  A(34,34) 

DIMENSION  FX(33,3  3),FY(33,33),GX<33,33),GY(33,33),UT(33,33),VT(3  3, 
333),  UAV(34,34),  VAV(34,34),  VORT ( 33  ,  3 3 )  , UAV E ( 34 , 34 ) ,  VAVE(34,34) 
DIMENSION  PS  IX  (35,33),  VORTX  (33,33) 
DIMENSI.L-N  VORT1  (33,33),  V0RT2  (33,33),  VORT3  (33,33) 

BETA  TERM  CALCULATION,  A/BETA  CONSISTENT  WITH  NORPAX  MODEL 

PI=2. 1415926 

OMEGA=2. *PI/86400. 

RAC=6. 375*10. **8 

BETA=(2. 0*OMEGA/PAD) *6  5./9  0. 

DELTAX=3.*10.**7 

WL=WICTH    OF    WESTERN    BOUNDARY;       WL 1=DEL TA X/2 ;        UNSTABLE    CASE 
WL1=DELT AX/2. 

ALFA  =  8ETA*(((SQf-T(3.)/(2.*PI))*WLl)**3) 
AMIN=ALFA 

IF     (KK.EG.O)     WRITE(6,28)     ALFA 
28    FORMAT (//, IX, 'ALFA    =»,E12.4) 
OX=DY    FDR    THIS    MODEL 
DELTAY=3.*10.**7 

CUP=.5 

DC  30  J=2,JMM1 
DO  30  1=2, IMM1 
30  Fid, J)  =  -(BETA)*(CUP*(PSI(I«-i,J)-PSHI,J))+(l.-CUP)*(PSI(I,J)- 
aPSH  1-1,  J)  )  )*DELTAX 

STRESS  CURL  TERM  CALCULATION  ********<<**************************<** 
THE  AMOUNT  OF  POSITIVE  VPRT1CITY  INTO  THE  OCEAN  FROM  THE  ATMOSPHERE 
IS  NEGATIVE  (  ANTI-CYCLONIC  ).  THIS  MAKES  THE  OCEAN  SPIN  A/C  LIKE  A 
HIGH, PRODUCING  A  POSITIVE  ST RE AMFUNC T I  ON ,  PS  I . 

F=l.O 

DEPTH=2. 0*10.0**5 
FJ*M1=JMM1 
B=FJMM1*DELTAX 
DO  40  J=2,JMM1 
Y=PI*(J-1)*(1.0/B)*OELTAX 

STRESS=  U2.*F*PI)  /  (DEPTH*3JI  *  SIN(2.*Y) 
DO  40  1=2, IMM1 
40  Fl  ( I ,J  )  =  F1(  I , J)-(STRESS  *DFLTAX  **2) 


DETERMI  JATIfN  OF  PSI  FIELD  BY  ANALYTICAL  MEANS. 
ANALYTIC  PSI  FIELD  WHICH  OCCURS  AT  ° OW  J=9. 


PRINT  OUT  MAX 


KK  =  KK«-1 

IF(KK.GE .2 )  GO  TO  43 

PSIANL  =-R  *  XIX)  *  f 


SETA**(-U)  *  CURL  (TAU) 


ANALYTICAL  PSI  FIELD 


P  =  DOMAIN  WIDTH 

FIMM1=  1.111 

R=F I MM1*( DELTAX) 

BETA  =  JF/DY  (AS  AB'TVE) 

CURL  (TAU)  =  K  COMPONENT  OF  WIND  STRESS  CURL  =  STRESS  (AS  ABOVE) 

XIX)  =  XX.   SEE  MUNK  (1951.  EQ.  20): 

XX=HKK*(EXP(-1/2*HK*XI )*COS( ( (SQRT  3  )/ 2 ) *HK* X+ ( SORT  3)/(2*HK*R)  - 
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42 
43 


PI/ 

HK  = 
HKK 
Y=P 
ST* 
J  =  9 
FOP 

del 

DEL 

DO 

X=< 

xx  = 

1*P  ) 

PSI 

nor 

PSI 

IP 


6)+l-( 1/HK*R) *(HK*X-EXP(-HK(P-X ) ) -1 ) ,     WHERE 

< BETA/ALFA)**. 33333 333 

=  2./SJRT     (3.)     -    SORK  ?.  )/ (HK*R) 

I*y.*(  l  ./B)*QELTAX 

ESS=-((2.*PI)/(DCPTH*R))*SIN<2.*Y) 


awRi 

44    FOR 

5)2,/ 

DEL 


EXPA  iU 
TAX  NOW 
TAX =3.* 
4  2  1=1, 
1-1  )*Ct 
-HKK* (L 
-PI/fj.) 
AMU  )  = 
MALIZfc 
ANL<  I  )  = 
( IMOIMK 
TE(6,44 
MAT  (  •  0' 
/ ilX, 33 
TAX=3.* 


ED  ANAL 

EQUAL 
10.**6 
100 
LTAX 
XP{ r.5* 
+  1.-U. 
-R*XX*6 
PSI  AM. 
PS  I  AM.  ( 
,20)  .EQ 
)  (PSIA 
,  ' THE  M 
F4.2,// 
10.**7 


YTICAL 
30  KM; 


SCLUTICN,  DFLTAX  DECREASED  BY  ORDER  OF  10 
WILL  EXAMINE  WESTERN  1/3  OF  OCEAN 


HK*X) )*CGS( ( (SQPT(3.) )/2. )*HK*X+( SORT ( 3. ) )/(2.*HK 
/(HK*R) )*(HK*X-EXP(-HK*(R-X) 1-1. ) 

ETA**(-1 ) *STPESS 

FIELD  3Y  DIVIDING  BY  PSIANL  MAX 

I  )/  (  L0.**8) 

.0)  .OR.  (KK.EQ.  1) ) 

NL(  I  )  ,1  =  1  ,99) 

AX    ANALYTIC    PSI     FIELD     IS    FOR    SOW    J=  9  :  '  ,  /  ,  I  X  ,  3  3F  4 . 

,  1X.33F4.2,// ) 


CALCULATION    OF    FRICTION    TERM    USING    PS1M1     (FORWARD    TIME    SCHEME)     * 

FRICTION  TERMS  ARE  CALCULATED  USING  A  FORWARD  TIME  SCHEIE.  THIS 
IS  DONE  TO  PRESERVE  LINEAR  COMPUTATIONAL  STABILITY  .IF  THE  LFAP 
FRCG    SCHEME    WAS    UTILIZED     ,     THE     FRICTION    TERM    WOULD    PE    UNSTABLE. 

DEFINF     INITIAL    U(32,32l     AND    V(32,32)     AS    FUNCTIONS    OF     PSIM1<33,33 
)  REFERENCE    SHOULD    BE    MADE    TO    SKETCH    OF    GRID     IN    NOTES:     ****** 

UU+1,J+1>  AND  V(H-l,J  +  l)  TRANSFORM  U,V(32,32)  INTERNAL  GRID  ** 
SYSTEM     INTO       EXTERNAL    U,V( 34,34)     GRID    SYSTEM    ******************* 
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50 


UO    50    1=1, IMM1 

DO    50    J=1,JMM1 

AA=    PS  IMK  1  +  1,  J+l) 

BB=    PS  IMK  I  ,  J+l)     + 

CC=    PSIMK  1  +  1,  JUI 

DD=    PS  IMK  Ifl,  J)     + 

EE     =    2.0*DELTAX 

V(I  +  1,J+1)     =     UAA/EE)     -(BB/EE)) 

UCI  +  liJ*-U     =     UDD/EE)     -(CC/EE)) 


+    PSIMKI  +  LtJ) 

PSIM1 ( I ,J) 

+  PS  IMK  I  ,  J+l) 

PSIM1 ( I ,J) 
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DEFINE    EXTERNAL    BOUNDARY    VALUES    OF     U    AND    V    ********************* 

DO    6  0     1  =  2.  IM 

U(  I  ,  1)     =     -U(  I ,2) 

V(  I  ,  1)     =    -V< I ,2) 

U(  I, JMPl  )=-U(  I ,JM) 

V(  I  ,  JMPl  )=-V( I,JM) 

DO     7C     J=2,JM 

U(  1, J)=-U(2, J  ) 

V(  1,  J) =-V(2, J) 

U<  IMP1 ,  J ) =-U( IV, J) 

V(  IMP1, J )  =  -V(  IV, J) 

U(  1,  1  I  =U(2,2  > 

V(  1,  1)=V(2,2J 

U(  34,34) =U( 33, 33  ) 

VI  34 ,34)  =V(33,33) 

U(  1, 34)=U( 2,33) 

V(  1,34>  =  V<  2, 33) 

LM'J4,1)=U133,2) 

VI 34,1)=V< 33,2) 

IF(KK.GE.O)     GO    TO 


3  ICO 


CALCULATION    OF     THE    FIELD    VARIABLE    OF    COEFFICIENTS    OF    EDDY    VISCOSITY 
•A,«     BASED    ON    ENSTROPHY    CASCADE 

V0PT(I,J)     IS    THE    RELATIVE    VOP.TICITY     AT     EACH    GRID    POINT. 
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1  101    DO    110f>     1=1, 

do  not)  j=i, 

1  lOt.    UAV(  I,  J)  =(U( 

on  mo   i  =  i, 

00    1110    J=l. 

1  110    VAV(  I, J) =( VI 

;  V0RT1CITY(I, 

00    1US     1  =  1, 

00    1115    J=l, 

1115    VJRT(I  ,J)=  (V 

3/0ELTAY 


IM 

JMI 

I 

I  Ml 

JM 

I 

J  ) 

IM 

JM 

AV(H-l.J)    -    VAV(  If  J))/DELTAX    -    (UAVU.J+1)    -    UAV(I.J)) 


IP1 

J) *U( I d,J) )/2. 
IP1 

I 

J) «-V( I ,J+1) )/2. 
=  i>V/DX    -    DU/OY 


IF  (K.GE.O  )  GO  TO  1120 


EXAMININ 
"NUMER IC 


C  VARIOUS  METHOOS  OF  CJMPUTING  VORTICITY  FROM  K.  MIYOKUOA 
AL  WEATHER  PREDICTION,  COMPUTATIONAL  METHOOS,"  1962. 


SETTING  UP  AN  EQUIVALENT  35X35  PSI  FIELU:  PSIX(I,J) 


121 


122 


123 


DO     112  1 
DO     1121 
PSIXd  ,  J 
DO     1122 
PSIX  (1,  J 
PSIX  (35, 
DO    1123 
P  S  I  X  (  I  ,  1 
PSIXd  ,3 
PSIXd,  1 
PSIX  (1,3 
PSIX(35, 
PSIX(35, 


1=2,34 
J=2, 34 

)=  psi(i-i,j-n 

J=2,34 
)=PSIX(3, J ) 
J) =PSIX( 33, J  ) 
1=2,34 
) =PSIX( 1,3) 
5)=PSIX( I ,33) 
)=PSIX(3,3) 
5)=PSIX( 3,33) 
1)=     PSIX(33, 3) 
35  )  =  PSIX(33, 33) 


ORIGINAL     EQUIVALENT     •X'    METHOD,    USES    5    GRIDPCINTS: 

DO  1124  1=2,34 
DO    1124    J=2,34 

124  V0RTX<I-1,J-1)=.5*IPSIX(I-1,J«-1)«-PSIX(I  +  1,J«-1)«-PSIX(I«-1,J-1)«- 
£PSIX(I-1,J-1)-4.*(PSIX(I,J)))/(DELTAX**2) 

MIYAKODA    METHOD    !,     'MINIMUM    GRIDPOINT'     SCHEME,     VORT+=VORTl,     USES 
5    GP.IDPDINTS 

DO  112  5  1=2,34 
DO    1125    J=2,34 

125  V0RT1(I-1,J-1)=    ,5*(PSIXU-1,J)«-PS.  dI,J«-l)*-PSlX(I«-l,J)*PSIX(I-l,J 
£)-4.*(P3IXU,J))  >/(DELTAX**2) 

MIYAKOCA    METHOD    II,     "0 *I FNTAT I  UN AL    MINIMUM    ERROR,"    VORT 2 = ( 2* VORT 1* 
l*V0RTX)/3;     USES    9    GP.IDPDINTS    TO    SOLVE     FOR    VORTICITY 

DO  1126  1=1,33 
DC  1126  J=l,33 
V0RT2( I, J) =(2.*V0RT1 ( I , J ) +VORTX ( I , J ) )/ 3. 

MIYAKOCA    METHOD    III,     "KNIGHTING    L    CGURA"     SCHEME:     < 2*VQRT 1-VCRTX ) /3 


1126 


DO    1127 
DC    1127 

1127  V0RT3( I , 

1128  DO    1129 
DO    1129 

1  129    VORTd  ,  J 
1120    CONTINUE 

DE TERM  II 
;  BASED    CN 

DO    1130 
DO    1130 
1130    DELVCkd 


1=1,33 
J=l,33 
J)=  (2.*V0PT1<  I,  J)  -VORTX(I.J)) 

1=1,33 
J=l,33 
)=V0RT3( I, J) 

ING  A  VARIA3LE  LATERAL  COEFFICIENT  OF  EDDY  V  I SCOS I T Y  ,  • A*  , 
ENSTROPHY  CASCADE 

1  =  2,  IM 
J=2, JM 
,  J)  =  SQRT(  ABS(  (  (  (VORT(  I  ,  J  )«-VORT  d , J-l ) ) /2 .-( VORT ( I-l  ,  JJ+VOR 
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o1T(  I-l.J-1)  )/2.)/DELTAX)**2M  ((  V0RTI1  ,J)+VORT  ( I  -1  ,  J  )  J  /2-  • 
«1 ) *VORT<  1-1,  J-l)  )/2.  )/DELTAX)**2H 

DVCRNX=2.*(10.**(-6)  )/DELTAX 
DO  1135  1=2, IM 
DO  1135  J=2,JM 


(VORT <  I,  J- 


1135    A(I,J)=    AMINM  l.+999.*DFLV0R(  I  ,  JJ/DVORMX) 

IFiKK.GE.O J    GO    TQ    1160 
:  BOUNCARY    VALUES    FOR    ENSTPQPHIC    A     (34,34): 

I  IN    ALL     FOIJR    DIRECTIONS 

DO    1150    1=2, IM 

SLCPE=(A(I  ,2)-A(  1,3)  )/DELTAY 

All,  1)  =  <  SLGPE*DELTAY)  +  A(  1,2) 

SLCPE=(A( I , JM)-A( I,JMM1) )/DELTAY 
11 50    A(I,JMD1)=(SL0PE*DELTAY)+A(I,JM) 

DO    1155    J=2,JM 

SLCPE=(A(2,J)-A(3,J) l/DELTAX 

All, J)-( SLGPE*DEITAX)  +  A<2,  Jl 

SLC°E=(A(If.J)-A(  I  MM  1, J)  )/DELTAX 
1155    A{  IMP1,.J  )=  (SLOPC*DELTAX)  *A(  IM,  J) 


LINEAR     EXTRAPOLATIONS 


1160  CONTINUE 
METHOD  2 
(  34,  34)  : 


FOR    DETERMINING    BOUNDARY    CONDITIONS    FOP     ENSTROPHIC     A, 
EXTENSION    OUTWARD    OF     EXISTING     INTERIOR    BOUNDARIES. 


1165 


1170 


DO  1165  J=2,JM 
A(  I  ,  J)  =t   (2,  J) 
A(IMP1,J)=A(IM,J) 
DO  1170  1=2, IM 
A(  I,  1) =A( I  ,2) 
A(  I  ,  JMP1 )=A(  I ,JM) 


FX=D/DX( A 
FY=D/DX( A 


CU/DX) 
CV/DX) 


0/DY( A 
D/DY{ A 


DU/DY) 
DV/DY) 


DO  1205  1  =  2, IM 

DO  1205  J=2,JM 

FX(I,J)  =  ((MI,JH-A(I«-1,J))/2.)*(U(H-1,J)-U(I,J))  - 
2UA(I-1,J)+A<I,J))/2.)*(U(I,J)-U(I-1,J))  + 
31 (  All, J) +A ( I , J  +  l )  )/2.  >*(U( I , J+l  )-U<  I, J)  )  - 
MlAl  I,  J-1)+A(I,J)  )/2.)*(U(I,  J)-UU,J-1)I 

FY(I,J)=((A(I,J)+A(I+1,J))/2.)*(V(I+1,J)-V(I,J)I  - 
2((A(I-1,J)+A(I,J))/2.)*(V(I,J)-V(I-1,J))+ 
J(  (  A(  I, J)  *A<  I  , J  +  1J  )/2 .  )*{ V( f  ,J*1)-V(  I, J)  )  - 
MlAl  I,  J-ll+AIIiJI  )/2.)*(V(  I,  J)-V(I,J-1)  I 

205  CONTINUE 
100  CONTINUE 

CALCULATION  CF  THE  FIFLD  VARIABLE  OF  COEFFICIENTS 

•A,'  BASED  CN  KINETIC  ENERGY  CASCADE 


OF  EDDY  VISCOSITY 


ENERGY  CASCADE,  DEFORMATION  ALONG  A  STR 
Y.  (DEFORMATION  NOPMAL  TO  A  STREAMLINE, 
LUTE  VALUE  OF  DEFORMATION,  DEF(I,J)  =  S 


FOR  KINETIC  E 
DV/DX  +  DU/DY 
DV/DY.   ABSC 
DEFT**2) . 

3101  DO  3105  1=1, IM 

DO  3105  J=1,JMP1 

VAVE(I,J)=(V(I,J)*V(I+l,J))/2. 
3  105  UAV(I,J)=(U(I,J)*U(I*l,J))/2. 

DO  3110  I=1,IMP1 

DO  3110  J=i,JM 

UAVE(I,J)=(U(I,J)+U(  I, J+l)  )/2. 
3110  VAVUtJ) -(V(!  ,J)+V(I  ,J  +  1)  )/2. 


A  STREAMLINE,  DEFS 

kltr,  DEFT  =  DU/DX 

SQRT(DEFS**2  * 


VAV  LSEO  IN  DV/DX,  UAV 
DAVE  JSED  IN  DU/DX 

DO  3  115  1  =  1,  I  M 
DO  3115  J=1,JM 


USED  IN  DU/DY,  VAVE  USED  IN  DV/DY, 
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DFFS( I  ,  J  ) 
<JAY 

DEFT(I,J  ) 
iJElTAY 
3115    DEF<  I, J)  = 


SCAL 
PSIM 
OEFM 

IF(K 
IF  (M 
GO  T 


E=10. 
AX  =  2. 
AX=2. 

K.ET. 
OD(Ki 
C  312 


=  (VAV(  1+1  ,  J  )-VAV(  I,  J)  » /DEL  TAX  MUAVI  I  ,  J*  1  I  - 
=  (UAVEl  1*1, JJ-UAVE(  I, J)  )/DELTAX-(VAVE( I  ,J+1 
SCRT(CEFS( I  ,J  >**2+UEFT( I , J) *«2) 

3*SCALE 
*PSIMAX/(I)ELTAX**2I 

1)  CO  TC  3121 

20)  .EU.O)  GO  TO  3121 
6 


UAV( I , J) l/DELT 
)-VAVE< I, Jl )/D 


3121 


DO  3 
JS=J 
DO  3 
RE5I 
WRIT 
FORK 
<i33(l 
3126  CONT 

AM  IN 
DO  3 
DO    3 


PRINT    OUT    ThE    DEFORMATION    FIELD 
=  lf  JM 


3122 
3125 


122  J 
MP1-J 
122     I 

C(  I  ,  J 
E<6,  3 
AT{  «0 
X,33F 

INUE 

=  ALFA 
130  I 
130    J 


=  1,  IM 

)=DEF(IfJS)/    DEFMAX 

125)     ((RESID    ( I , J  )  ,  1  =  1  ,  33) , J=l,33) 

■  ,'P.ESID    FI ELD=DEFORMATION/DEFORMATION    MAX* 
4.3,//) ) 


=  1,  IM 

=  1,  JM 


,/, 


3130    AU,J)=AMIN*(  l.*999.*DEF<  I,  J  I /DEFMAX  J 

IF(KK.GE.O)    GO    TO    3210 


:  FX=D 

;  FY=D 

DO    3 

DO    3 

FX(I 

2((  A( 

3(  (A( 

4((A{ 

FY(I 

2(  (  A( 

3(  (A( 

4(  (A( 

3205    CCNT 

GO  T 
3210    CCNT 

FXI! 
2((A( 
3(  (A( 
M(A{ 

FY(I 

2(  (A( 

3(  (  A( 

4(  (A( 

3215    CONT 

CALC 
INTE 

CURL 

DO  1 
DO  1 
CURL 

a-(  (F 

Fl  (I 
1210    CC\T 


/DX(  A 
/DX(  A 

205  I 
205  J 
,J  )=  ( 
1-1,  J 
li  J)  ♦ 
I,  J-i 

,J)  =  ( 
I-lt  J 
I*  J)  + 
It  J-l 

IKUE 

0  321 
INUE., 
,J)=(( 

1  -  i  t  J 
I  ,  J)  ♦ 
1,  J-l 


DU/DX)     ♦    D/DY(A    DU/DY) 
CV/DX)     +    D/DY(A    DV/DY) 

=  2,  IM 
=  2,  JM 

(MIiJ-l)«-AU.J))/2.)*(U(H-l,J)-U(I,jn     - 
-1)+A(I-1,J)     )/2.»*<U(I,J)-U(I-l,Jll«- 
A(I-l,J))/2.     )*(U(  I  ,J*1)-U(  I.  J)  )     - 
MA(I-l,J-l))/2.)*(U(ItJ)-U(I,J-l)) 

<AU,J-1)*A(I,J))/2.)*(V(I-H,J)-V(I,J))     - 
-1  J  +A( I -1 , J  »     )/2.)*<V(I.J)-V<  1-1,  J)  )* 
A(I-l,J))/2.     )*(V(I  ,  J  +  li-V(  I, J)  )     - 
)+A(I-l,J-l))/2.)*(V(I,J)-V<I,J-l)) 


,JK< 
I  - 1  ,\J 
I,  jW 
It  J-l 

INUE 

ULATE 

RIOR  ( 

F    = 

210  I 
210  J 
F=  {(  F 
X(  1-1 
-ltJ- 
INUE 


TA(I,J-l)+A(l.J))/2.)*((DEFT(l,j)+DEFT(I,J- 
-1) +A ( I -1  ,  J  )  j/2.  )* (  ( DEFT( 1-1 ,J )>OEFT ( I-l, J 
A(I-l,J))/2.  )*<(DEfS<I,JH-DEFS(J-l,J))/2.) 
)*A(  1-1, J-l ) )/2. ) *l  (DEFS(I,J)t-L)EFS(I-l,J-l) 

<A(I,J-1)+A(I,J))/2.)*((DEFS(I,J)*-QEFS(I,J- 
-1)*A(I-1,J)  )/2.)*((DEFS(I-l,J>*DEFS<I-l,J 
AU-l,J))/2.  )*(  (DEFT  (  I  ,  J  J+DEFT  (  1-1  ,J  M/2.J 
)+A(I-l,J-l))/2.)*({DEFT(I,J)*aEFT(I-l,J-l) 

THE    CURL    OF    THE    FRICTION    FORCE,     DEFINED    ON 
31,31) 
CFY/DX    -    DFX/DY 

=  3,  IM 

=  3,  JM 

Y(  I  ,  J)  +  FY(  I  ,J-l)  )/2.    -     (FY(  1-1,  J)«-F  Y{  1-1,  J- 

,J)*FX(I,J))/2.-(FX(I-l,J-lH-FX{I,J-l))/2.) 

1)=CURLF     *F1( 1-1, J-l) 


Do  321$  1-2, IM 
Do  3215  J-2,JM 

1)  1/2. )- 
-1)  J/2.)* 

)/2.))*DSLTAX 

1)1/2. )- 
-1)  )/2.)  + 

,/2-O^DELTAX 

THE  GRID 


1  )  )/2. I/DELTAX 
/DELTAY 


DO  150  I=1,IMM1 
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00  150  J=l,JMMl 

AA=PSI  (  I  H,JU  >♦  PSI  (  I  «-l  ,  J  ) 

BB=PSI  <  I  ,  J«-1)+PSI  U  ,  J) 

CC  =  PSI  (  I  4-1 ,  J  +  l  )«-PSI  <  I.  J*l) 

D0=PSI(I  «-l,  J)*PSH  I,  J) 

EE=2.0*DEL  TAX 

UT(I ,J)=0. 

VT( I  ,J)  =  0. 

V(  !, J)  =  ( AA/tE J-( BB/EE) 
150  U(I(J)>( UD/EE)-(CC/EE)) 
C 

C      EASTWARD  MASS  FLUX 
C 

00  4  00  I =1  ,  IKM1 

FX(I ,1)=U<  I,  1) 

DO  300  J=2,JMM1 
3  00  FX(l,J)=0.5MU(I,J)«-U(ItJ-l)J 
400  FX( I  ,JM)=U( ! »JMM1) 
C 

C      NORTHWARD  MASS  FLUX 
C 

DO  6C0  J=1,JMM1 

FY<  1, J)=V(  1, J) 

DO  500  I=2,IMM1 
500  FY{ I  ,J)=0.5*(V(I , J)*  V( 1-1 1 J ) J 
600  FYUN,  J)  =V(  IMM1, J) 
C 

C      EASThAR}  ACVECTION  OF  MOMENTUM  V=(FX,FY) 
C 

FACT=0.125 

DO  800  !=2,!MM1 

IM1=I-1 

DO  700  J=1,JM 
700  GX(I  t  J)  =  FX(  I  ,  J)«-FX(IM1,J) 

DO  800  J=1,JMM1 

FLLX=FACT*<GX(  I  ,  J ) *GX ( I , J+ 1 )) 

UFLLX=CJ(  I  ,  Jl+Ul  I  Ml,  J)  )*FLUX 

UT(  I  ,  J  )=UT  (  I  ,  J  J+-UFLUX 

UTIIM1,  Jl^UTdfl,  Jl-UFLUX 

VFLUX=IV( I fJ)+V( I  Ml, J) )*FLUX 

VT(  I  ,J»=VTi  I,  J)«-VFLUX 
800  VTUM1,  J)  =  VT(  IM1  ,  J)-VFLUX 
C 

C      NOPTHrtARC  ADVECTION  OF  MOMENTUM  ( D/ DY( VU ) *D/ DY ( V V ) ) 
C 

DO  1C00  J=2,JMV1 

JM1=J-1 

DO  9C0  1=1, IM 
900  GY(I,J)=FYII  ,  JM1  J4-FY(I  ,J  ) 

DO  1000  I=1,IMM1 

IP1=I+1 

FLLX=(GY( I ,J)+GY( I  PI , J ) ) *FACT 

UFLLX=(U( I  ,  J)*U(  I , JM1)  )*FLUX 

UT(1  ,J)=bT<  I  ,  J)+UFLUX 

UT(I,JM1)  =  UT(I,  JMD-UFLUX 

VFLUX=(V(I , J)+V(I, JM1) )*FLUX 

VT(  I  ,  J  )  =  VT  (  I  ,  JH-VFLUX 
1C00  VT(  I  ,JM1  )  =  VT<  I  f  JMD-VFLUX 
C 

C      CURL  OF  AUVECTICN  TERMS 
C 

DO  190  I=2,IMM1 

DO  190  J=2,JMM1 

AAA=(VT(  I  ,J) fVT[ I . J-l)  )/2. 

BBB=(VT(  I-l.Jl*VT{  I-l.J-1)  )/2. 

CCC=(UT(  I  , J)+UT{ I  - 1 , J) )/2. 

CDD=<UT( I  ,  J-1)+UT(  I-1.J-1M/2. 

DELVX=  AAA  -  BB6 

DELUY=  CCC  -  DDO 
190  F  1  (  I  ,  J ) =  F1(I|J)+   DELVX  -  UELUY 
C 

C      SET  LP  FCR  SOLVER  (UPSIDT=0  ON  BOUNDARIES 
C 

DO  2CC0  J=2,JMK1 
DO  2000  1  =  2, ^Ml 
2000  DPSICTt I , J)=-Fl(  I,  J) 

RETURN 
END 
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SURPCUTINE     SOLVE R(  B,  M  ,  N,  OF LTAX , DELT AY) 
FAULKNER        DECK    QQ 

RCLANO    SWEET'S    PCISSON    SOLVER    FOR    A    RECTANGLE 
REVISED    AUG    I r,  74     FOR    PROF     ROBERT    HANEY 

THIS     SUiPPOTINE    SULVES    THE     EQUATION    LAPLACIAN(Q)     =    F(I,J),     AND    RETURNS 

THE     SCLUTICrj     IN    R(  I  ,  J)  . 

B    =    eiM*-  l.N+1  )  ,     WHERE     N    =    A    POWER    OF    T  WC  . 

B    =    -F*0ELTAY**2     (UNITS    OF    Q)     CM     IMTEPIP"    POINTS. 

B    =     SPECIFIED    VALUES    OF    Q    ON    BOUNDARIES. 

LAPLAC  IAN(  r.)     =    5-POIMT     DIAMOND    APPRPX. 

DIMENSUN    B(33,33)f     P(33),     TWOCOS(33),     ^ECIP(32) 


NP1=N+1 


c 

c 

THIS    SECTICN    GENERATES    THE    ROOTS    CF    THE    POLYNDMIALS 

C 

c 

FOP    THE     REDUCTION    AND    SlLUtION 

LO    =    N/2 

TWCCOS(lO)        =    0. 

110 

L     =     LD/2 

TWUCPS(L)     =    SURH2.0    +    TWOCOS(LO)) 
LO    =    L 

120 

TWCCOSIN-L)     =    -TWOCOS(L) 

L     =     L    +    2*L0 

IF     (  (2*L/N)*(2*L0-3)  )     14  0,130,110 

130 

TWOCOS(L)     =     (Tr-JCCS(  L+LO)     ♦■    TWOCOS ( L-L 0 )) /TWOCOS ( LO ) 
GO    TO     1J0 

C 

140 

CCNT  INUE 

S    =     (DlLTAX/0ELTAY)**2 

LO    =    N/2 

JU    =    N 

IU    =    M-l 

MM2     =    M    -    2 

JUKI    =     JU    -     1 

P(")     =    0. 

RIGHT    HAND    SIDE     IS    MULTIPLIED    BY    DELTAY**2,     NOWAOD    THE    PROPER     MULTIPLI 
OF    THE    UPPER     AND    LOWER    BOUNDARY    DATA.       THE    SIDE    BOUNDARY    DATA     IS 
ENTERED     INTO    THE    CALCULATION    DURING     THE    SEDUCTION. 


200 


220 


300 


B(  1,  J)/S 
BIM+1, J)/S 


LI/N 


DO    2C0    J=2,JU 
B(2, J)     =       B(2, J>     • 
B(Mt  J)     =        BC,  J) 

START    PEDUCTICN 

ID    =     1 

MODE     =     2 

LI     =     2*L0 

IPHASE    =    2*M0DE    - 

JO    =    N/L I 

JH    =    N/( 2*LI  1 

JT    =     JD    +     JH 

JI     =    2*JD 

JO    =    JD*MODL     ♦    1 

DO     750    J=JC,JU,JI 

IPHASE    =  3    FOR  THE  INITIAL     STEP    OF    THE     REDUCTION 

=  t*    FOR  THE  REMAINING    STEPS    OF    THE     REDUCTION 

=  2    FUR  THE  FIRST     STEPS    OF    tH£     SOLUTION 

=  1     FCR  the  LAST    STEP    CF     THE     SOLUTION    BEFORE     EXITING 

GO    TO     (370,350,330,  310),  IPHASE 

ACCORDING  TO  IPHASE  SET  UP  PROPER  P I GH T  HAND  SIDE  (P  ARRAY)  AND 
AMOUNT  TO  BE  ADDED  ( J-TH  COLUMN  OF  B  AR^AY)  TO  THE  SOLUTION  OF 
A  CERTAIN  SET  OF  TRIDI AGONAL  SYSTEMS 
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310 

320 
330 
340 
350 
360 
370 
380 


400 


500 

600 
700 

750 

800 


840 
850 

900 


DO  3  20  I  =2,M 

PI  =  B( I  ,J  )  -  B(  I, 

8( I  ,  J)  =  B( I , J)  - 

P(  1-1)  =  B( I  ,J)  + 

GO  TO  40C 

DO  340  I  =2  ,M 

P{ 1-1)  =2.*B(  I  ,J) 

B(  I, J)  =   B(  I,  J*l) 

GO  TO  400 

00  360  I=2,M 

P(l-l)  =    B  {  I  ?  J  ) 

B< I  ,  J)  =(B(  I  ,  J)  - 

GO  TC  4C0 

DO  3B0  1=2, M 

P(  1-1)  =    B(  I, J) 

B(  I,  J)  =  0. 


J-JT)  -  B<  I  ,  J*JT) 


B(  I 

PI 


J-JH)  -  BUtJ+JH}  +  BU.J-JD)  ♦  B(I,J«-JD) 


♦  B(I ,J-l) 


+  8(1, J-JD) 
B{ I, J-JH)  - 


<■  B(  I  , J+JD) 
Bl  I, J  +  JH)  )/2. 


♦  B( I, J-l)   +  B( I , J+ll 


SOLVE  A  TPIPIAGONAL  SYSTEM  WHOSE  COEFFICIENT  MATRIX  HAS  THE  FORM 
TRIDIAG(-1,2  «■  S*(2  -  TwrCOS(D),  -1)  AND  RIGHT  HAND  SIDE  IN 
P  AFRAY.   SOLUTION  IS  BY  GAUSS  ELIMINATION. 

DO  7C0  L=LO,N,LI 

A  =  2.+  S*(2.-  TWOCOS(L) ) 

RECIP( IU)  =  -l./A 

P(  IU)  =  P(  IU)*S/A 

DO  503  IP=2,MM2 

I  =  P  -  IP 

RECIP(I)  =  -l./( A+RECIP(  Itll  ) 

P(I)  =  (S*P(I)  ♦  P(  1+1  )  )/( A+RECIP<  I  +  l)  ) 

P(l)  =  (S*P(1)  +  P(2   )  )/( A*RECIP<2 )  ) 

DO  6C0  1=2,  IU 

P(  I  )  =  P( I  )  -  RECIP( I )*P(I-1) 

CONT  INUE 

DO  750  1=2, M 

8(1  , J)  =  B(I  ,J)  +  P(  1-1) 

GO  TO  1900,850,800,800),  IPHASE 

LO  =  LO/2 

MODE  =«  1 


IF  (LO  .EQ. 
GO  TO  220 
IF(LO.GT.l) 
LO  =  2*L0 
IF  (LO  .LT. 
CONTINUE 
RETURN 
END 


1) 


GO   TO  220 

N)  GO  TO  220 
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